#include<math.h>
#include<iostream>
#include <direct.h>
#include<ctime>
using namespace std;

double* initCondition(double xl, double xr,int N,double pe,double Ma)
{
    double dx,x[N+1],u[N+1],p[N+1],rho[N+1],Ax[N+1];
    int i,j;
    dx = (xr - xl) / N;
    for(i=0;i<N+1;i++){
   	x[i] = xl + dx * i;
	Ax[i] = 0.5 + pow(x[i],2);
    }
    for(i=0;i<N+1;i++){
	u[i] = 2 ;
	p[i] = 2 ;
	rho[i] = 2 ;
    }

	u[0] = Ma * pow((1.4 / (1 + 0.2 * pow(Ma,2))),0.5);
	u[N] = 2 * u[N-1] - u[N-2];
	p[0] = 1 / pow((1 + 0.2 * pow(Ma,2)),(1.4 / 0.4));
	p[N] = pe;
	rho[0] = 1 / pow((1 + 0.2 * pow(Ma,2)),(1 / 0.4));
	rho[N] = rho[N-1] * pow((p[N] / p[N-1]),(1 / 1.4));
    
    double out[6][N+1];
    for(i=0;i<N+1;i++){
        out[0][i]=u[i];
        out[1][i]=p[i];
        out[2][i]=rho[i];
        out[3][i]=x[i];
        out[4][i]=Ax[i];
        out[5][0]=dx;
        }
    double* pos;
    pos=&out[0][0];

	return pos;
}

int Laval_precise(double xl, double xr, int N, double pe, double Ma)
{
    double dx,x[N+1],p[N+1],Ax[N+1],ma[N+1],Maa[N+1],A1=0.5;
    int i,j;
    dx = (xr - xl) / N;
    for(i=0;i<N+1;i++){
   	x[i] = xl + dx * i;
	Ax[i] = A1 + pow(x[i],2);
    }

    int pos=1;
    while(pos<N){
        for(i=0;i<N+1;i++){
        ma[i]=1;
        Maa[i]=1;
    }
        for(i=0;i<pos;i++){
            double err=1;
            Maa[i]=Maa[i]+1;
		while(err>0.0001){
			ma[i]=pow((1.2*pow((pow(Maa[i],2)*pow(Ax[i]/A1,2)),(0.4/2.4))-1),0.5)*pow(2/0.4,0.5);//sup
			err=fabs(ma[i]-Maa[i]);
			Maa[i]=ma[i];
        }
		p[i] = 1 / pow((1 + 0.2 * pow(Maa[i],2)),(1.4 / 0.4));
    }

    double M2,p2,p02,A2;//normal wave
    M2=pow(((1+0.2*pow(Maa[i-1],2))/(1.4*pow(Maa[i-1],2)-0.2)),0.5);
	p2=p[i-1]*(1+(2.8/2.4)*(pow(Maa[i-1],2)-1));
	p02=p2*pow((1+0.2*pow(M2,2)),(1.4 / 0.4));
	A2=A1/p02;

    for(i=pos;i<N+1;i++){
            double err=1;
            Maa[i]=Maa[i]-0.5;
		while(err>0.0001){
			ma[i]=(A2/Ax[i])*pow((1+0.2*pow(Maa[i],2))/1.2,(2.4/0.8));//sub
			err=fabs(ma[i]-Maa[i]);
			Maa[i]=ma[i];
        }
		p[i] = p02 / pow((1 + 0.2 * pow(Maa[i],2)),(1.4 / 0.4));
    }

    pos++;

    if(fabs(p[N]-pe)<0.001){
        break;
    }

    std::string prefix = "./Output";
    if (_access(prefix.c_str(), 0) == -1)
        _mkdir(prefix.c_str());    
    FILE* fpout;
    char str[20];
    sprintf(str, "./Output/precise.dat");
    char *path=str;
    fpout = fopen(path, "w+");
    fprintf(fpout, "variables=\"x/L\" \"Ma\"\n");
    fprintf(fpout, "ZONE T=\"ZONE 001\" i=%d F=POINT\n",N+1);
    for (int i = 0; i < N+1; i++){
        fprintf(fpout, "%15.8f\t%15.8f\t\n", x[i],Maa[i]);  //把数据写入文本
    }
    fclose(fpout);

    }

	return pos;
}

double Find_min(double *ary, int left, int right, int n)
{
    double array[n];
    for(int j=0;j<n;j++){
        array[j]=*((double*)ary + j);
    }

	double min=fabs(array[left]);
	for(int j=left;j<right;j++){
		if(fabs(array[j])<min){
            min=fabs(array[j]);}
    }
	return min;
}

double Find_max(double *ary, int left, int right, int N)
{
    double array[3][N+1];
    for(int j=0;j<3;j++){
        for(int i=0;i<N+1;i++)
            array[j][i]=*((double*)ary + j*(N+1)+i);
    }

	double max=fabs(array[0][left]);
    for(int j=0;j<3;j++){
	    for(int i=left;i<right;i++){
		    if(fabs(array[j][i])>max){
                max=fabs(array[j][i]);}
    }
    }
	return max;
}

double calTimeStep(double dx, double CFL, double* pos_u, double* pos_a,int len){
    double u[len],a[len],dt[len];
    for(int i=0;i<len;i++){
        u[i]=*((double*)pos_u+i);
        a[i]=*((double*)pos_a+i);
    }
    for(int i=0;i<len;i++){
	dt[i] = CFL * dx / (a[i] + fabs(u[i]));
    }
    double* pos;
    pos=&dt[0];
    double min;
    min=Find_min(pos,0,len,len);
	return min;
}

double* energy(double* pos_u, double* pos_p, double* pos_rho,int N)
{
    double u[N+1],p[N+1],rho[N+1];
    for(int i=0;i<N+1;i++){
        u[i]=*((double*)pos_u+i);
        p[i]=*((double*)pos_p+i);
        rho[i]=*((double*)pos_rho+i);
    }

    double out[N+1];
    for(int i=0;i<N+1;i++){
        out[i]=1 / 0.4 * p[i] / rho[i] + 0.5 * pow(u[i],2);
    }
    double* pos;
    pos=&out[0];
	return pos;
}

double* deEnergy(double* pos_U, int N)
{	
    double U[3][N+1],out[N+1];
    for(int i=0;i<N+1;i++){
        U[0][i]=*((double*)pos_U+i);
        U[1][i]=*((double*)pos_U+(N+1)+i);
        U[2][i]=*((double*)pos_U+2*(N+1)+i);
    }
    for(int i=0;i<N+1;i++){
    out[i]=0.4 * (U[2][i] - 0.5 * pow(U[1][i],2) / U[0][i]);
    }
    double* pos;
    pos=&out[0];
    return pos;
}

class Laval
{
public:
		double Cx,CFL,dx;
        int N;
		double* Couples(double* pos_u, double* pos_p, double* pos_rho, double* pos_x, double* pos_Ax, int N);
        double* deCouples(double* pos_U,int N);
        double* macCormack(double* pos_uprxAx, double dx, double CFL, double Cx, double errlimit, int iterlimit, int N);
};

	double* Laval::Couples(double* pos_u, double* pos_p, double* pos_rho, double* pos_x, double* pos_Ax, int N)//ok
    {

    double u[N+1],p[N+1],rho[N+1],x[N+1],Ax[N+1];
    for(int i=0;i<N+1;i++){
        u[i]=*((double*)pos_u+i);
        p[i]=*((double*)pos_p+i);
        rho[i]=*((double*)pos_rho+i);
        x[i]=*((double*)pos_x+i);
        Ax[i]=*((double*)pos_Ax+i);
    }
    double *posu,*posp,*posr,*posx,*posAx;
    posu=&u[0];
    posp=&p[0];
    posr=&rho[0];
    posx=&x[0];
    posAx=&Ax[0];

    double* pos_E;
    double E[N+1];
	pos_E = energy(posu, posp, posr,N);
    for(int i=0;i<N+1;i++){
        E[i]=*((double*)pos_E+i);
    }
    
    double U[3][N+1],F[3][N+1],H[3][N+1];
    for(int i=0;i<N+1;i++){
		U[0][i] = rho[i];
        U[1][i] = rho[i] * u[i];
        U[2][i] = rho[i] * E[i];
    }
    double *posU,*posde;
    posU=&U[0][0];
    posde=deEnergy(posU,N);
    double deEn[N+1];
    for(int i=0;i<N+1;i++){
            deEn[i]=*((double*)posde+i);
    }

    for(int i=0;i<N+1;i++){
        F[0][i] = U[1][i];
        F[1][i] = pow(U[1][i],2) / U[0][i] + deEn[i];
        F[2][i] = (U[2][i] + deEn[i]) * U[1][i] / U[0][i];
        H[0][i] = 2 * x[i] * U[1][i] / Ax[i];
        H[1][i] = 2 * x[i] * (pow(U[1][i],2) / U[0][i]) / Ax[i];
        H[2][i] = 2 * x[i] * ((U[2][i] + deEn[i]) * U[1][i] / U[0][i]) / Ax[i];
    }

    double out[9][N+1];
    for(int i=0;i<N+1;i++){
        out[0][i] = U[0][i];
        out[1][i] = U[1][i];
        out[2][i] = U[2][i];
        out[3][i] = F[0][i];
        out[4][i] = F[1][i];
        out[5][i] = F[2][i];
        out[6][i] = H[0][i];
        out[7][i] = H[1][i];
        out[8][i] = H[2][i];
    }
    double* pos_out;
    pos_out=&out[0][0];
    return pos_out;
    }

    double* Laval::deCouples(double* pos_U,int N)
    {
    double u[N+1],p[N+1],rho[N+1];
    for(int i=0;i<N+1;i++){
        rho[i]=*((double*)pos_U+0*(N+1)+i);
        u[i]=*((double*)pos_U+1*(N+1)+i)/rho[i];
        p[i]=0.4*(*((double*)pos_U+2*(N+1)+i)-0.5*rho[i]*pow(u[i],2));
    }
    double out[3][N+1];
    for(int i=0;i<N+1;i++){
            out[0][i]=u[i];
            out[1][i]=p[i];
            out[2][i]=rho[i];
    }
    double* pos;
    pos=&out[0][0];
	return pos;
    }

    double* Laval::macCormack(double* pos_uprxAx, double dx, double CFL, double Cx, double errlimit, int iterlimit, int N)
    {
        double u[N+1],p[N+1],rho[N+1],a[N+1],Ax[N+1],x[N+1];
        int i;
        for(i=0;i<N+1;i++){
        u[i]=*((double*)pos_uprxAx+0*(N+1)+i);
        p[i]=*((double*)pos_uprxAx+1*(N+1)+i);
        rho[i]=*((double*)pos_uprxAx+2*(N+1)+i);
        x[i]=*((double*)pos_uprxAx+3*(N+1)+i);
        Ax[i]=*((double*)pos_uprxAx+4*(N+1)+i);
        }

        double* pos_u,*pos_a,*pos_p,*pos_rho,*pos_xx,*pos_Axx;
            pos_u=&u[0];
            pos_a=&a[0];
            pos_Axx=&Ax[0];
            pos_xx=&x[0];
            pos_p=&p[0];
            pos_rho=&rho[0];
        int times=0;
        double err=1;
        
		do{
            for(int i=0;i<N+1;i++){
            u[i]=*((double*)pos_u+i);
            p[i]=*((double*)pos_p+i);
            rho[i]=*((double*)pos_rho+i);
            }
            for(int i=0;i<N+1;i++){
            a[i]=pow(1.4*p[i]/rho[i],0.5);
            }
            double dt;
			dt = calTimeStep(dx, CFL, pos_u, pos_a, N+1);
            double*pos_UFH;
			//Predictive Step
			pos_UFH = Couples(pos_u, pos_p, pos_rho, pos_xx, pos_Axx, N);
            double U[3][N+1],F[3][N+1],H[3][N+1];
            for(int i=0;i<N+1;i++){
                U[0][i]=*((double*)pos_UFH+0*(N+1)+i);
                U[1][i]=*((double*)pos_UFH+1*(N+1)+i);
                U[2][i]=*((double*)pos_UFH+2*(N+1)+i);
                F[0][i]=*((double*)pos_UFH+3*(N+1)+i);
                F[1][i]=*((double*)pos_UFH+4*(N+1)+i);
                F[2][i]=*((double*)pos_UFH+5*(N+1)+i);
                H[0][i]=*((double*)pos_UFH+6*(N+1)+i);
                H[1][i]=*((double*)pos_UFH+7*(N+1)+i);
                H[2][i]=*((double*)pos_UFH+8*(N+1)+i);
            }

            double dU_dt[3][N+1],S[3][N+1],Us[3][N+1];
            for(int j=0;j<3;j++){
                dU_dt[j][0]=0;
                dU_dt[j][N]=0;
                S[j][0]=0;
                S[j][N]=0;
            }
            for(int j=0;j<3;j++){
                for(int i=1;i<N;i++){
                    dU_dt[j][i] = -((F[j][i+1] - F[j][i]) / dx + H[j][i]);
                    S[j][i] = Cx * fabs(p[i+1] - 2 * p[i] + p[i-1]) / \
                        (p[i+1] + 2 * p[i] + p[i-1]) * (U[j][i+1] - 2 * U[j][i] + U[j][i-1]);	// Artificial Viscosity Add
                }
            }
            for(int j=0;j<3;j++){
                for(int i=0;i<N+1;i++){
                    Us[j][i] = U[j][i] + dU_dt[j][i] * dt + S[j][i];
                }
            }

			double* pos_Us,*pos_Usout;
            pos_Us=&Us[0][0];
            double us[N+1],ps[N+1],rhos[N+1];
            pos_Usout= deCouples(pos_Us,N);
            for(int i=0;i<N+1;i++){
                us[i]=*((double*)pos_Usout+0*(N+1)+i);
                ps[i]=*((double*)pos_Usout+1*(N+1)+i);
                rhos[i]=*((double*)pos_Usout+2*(N+1)+i);
            }
			us[N] = 2 * us[N-1] - us[N-2];
			rhos[N] = rhos[N-1] * pow((ps[N] / ps[N-1]),(1 / 1.4));

			// Correction Step
			double* pos_ps,*pos_rhos,*pos_us,*pos_UFHs;
            pos_us=&us[0];
            pos_ps=&ps[0];
            pos_rhos=&rhos[0];
            
			//Predictive Step
			pos_UFHs = Couples(pos_us, pos_ps, pos_rhos, pos_xx, pos_Axx, N);
            for(int i=0;i<N+1;i++){
                Us[0][i]=*((double*)pos_UFHs+0*(N+1)+i);
                Us[1][i]=*((double*)pos_UFHs+1*(N+1)+i);
                Us[2][i]=*((double*)pos_UFHs+2*(N+1)+i);
                F[0][i]=*((double*)pos_UFHs+3*(N+1)+i);
                F[1][i]=*((double*)pos_UFHs+4*(N+1)+i);
                F[2][i]=*((double*)pos_UFHs+5*(N+1)+i);
                H[0][i]=*((double*)pos_UFHs+6*(N+1)+i);
                H[1][i]=*((double*)pos_UFHs+7*(N+1)+i);
                H[2][i]=*((double*)pos_UFHs+8*(N+1)+i);
            }
			double dUs_dt[3][N+1],dU_dtav[3][N+1],Ss[3][N+1];
            for(int j=0;j<3;j++){
                for(int i=0;i<N+1;i++){
                    dUs_dt[j][i] =0;
                    dU_dtav[j][i] = 0;
                    Ss[j][i] = 0;
                }
            }
            for(int j=0;j<3;j++){
                for(int i=1;i<N;i++){
                    dUs_dt[j][i] = -((F[j][i] - F[j][i-1]) / dx + H[j][i]);
                    dU_dtav[j][i] = 0.5 * (dUs_dt[j][i] + dU_dt[j][i]);
                    Ss[j][i] = Cx * fabs(ps[i+1] - 2 * ps[i] + ps[i-1]) / \
                    (ps[i+1] + 2 * ps[i] + ps[i-1]) * (Us[j][i+1] - 2 * Us[j][i] + Us[j][i-1]);		// Add artificial viscosity again
                }
            }
            for(int j=0;j<3;j++){
                dU_dtav[j][0] = 0.5 * (dUs_dt[j][0] + dU_dt[j][0]);
                dU_dtav[j][N] = 0.5 * (dUs_dt[j][N] + dU_dt[j][N]);
            }

            double Ut[3][N+1];
            for(int j=0;j<3;j++){
                for(int i=0;i<N+1;i++){
                    Ut[j][i] = U[j][i] + dU_dtav[j][i] * dt + Ss[j][i];
                }
            }

			double* pos_Ut,*pos_Utout;
            pos_Ut=&Ut[0][0];
            double ut[N+1],pt[N+1],rhot[N+1];
            pos_Utout= deCouples(pos_Ut,N);
            for(int i=0;i<N+1;i++){
                ut[i]=*((double*)pos_Utout+0*(N+1)+i);
                pt[i]=*((double*)pos_Utout+1*(N+1)+i);
                rhot[i]=*((double*)pos_Utout+2*(N+1)+i);
            }
            
			ut[N] = 2 * ut[N-1] - ut[N-2];
			rhot[N] = rhot[N-1] * pow((pt[N] / pt[N-1]),(1 / 1.4));
            for(int i=0;i<N+1;i++){
			u[i]= ut[i];
            p[i]= pt[i];
            rho[i]= rhot[i];
            }
            pos_u=&u[0];
            pos_p=&p[0];
            pos_rho=&rho[0];

            double* pos_dU_dtav;
            pos_dU_dtav=&dU_dtav[0][0];
            err=Find_max(pos_dU_dtav, 1, N, N);
            times++;
			if (times==iterlimit){
				break;
			}
		}while(err>errlimit);
    double out[3][N+1];
    for(int i=0;i<N+1;i++){
            out[0][i]=u[i];
            out[1][i]=p[i];
            out[2][i]=rho[i];
    }
    double* pos;
    pos=&out[0][0];
	return pos;
    }

int main()
{   
    clock_t startTime,endTime;
	startTime = clock();
    cout << "Program is running\n" <<endl;

    //read frome Input.txt
	FILE *fp = NULL;
    char *file = (char*)"input.txt";
    char line[100];
    if( (0 != access(file,R_OK|F_OK)) || (NULL==(fp=fopen(file,"r"))) )
    {
        printf("open %s failed\n",file);
        return -1;
    }
    while( fgets(line, 256, fp) != NULL ){}
    if(fp!=NULL)
    {
        fclose(fp);
    }
	const char *d = ",";
    char *pp,*inti;
	int kk=0,inp[2];
	double output[9];
    pp = strtok(line,d);
    while(pp)
    {
		char str[10];
		string ing=pp;
		int ii;
    for( ii=0;ii<ing.length();ii++)
        str[ii] = ing[ii];
    str[ii] = '\0';
		double dd;
		output[kk]=atof(str);
        pp=strtok(NULL,d);
		kk++;
    }
    double Ma=output[0],xl=output[1],xr=output[2],pe=output[3],Cx=output[4],CFL=output[5],errlimit=output[6];
    int i=0,iterlimit=output[7],N=output[8];
    cout << "Input:\nMa=" << Ma << "\tx_left="<< xl<< "\tx_right="<<xr<< "\tpe="<<pe<< "\tCx="<<Cx\
        << "  \tCFL="<<CFL<< "   \tErrlimit="<<errlimit<<"   \tIterlimit="<<iterlimit<<"   \tN="<<N<<'\n' <<endl;

    double *init;
    init=initCondition(xl, xr, N, pe, Ma);
    double dx,x[N+1],u[N+1],p[N+1],rho[N+1],Ax[N+1];
    for(int i=0;i<N+1;i++){
        u[i]=*((double*)init+0*(N+1)+i);
        p[i]=*((double*)init+1*(N+1)+i);
        rho[i]=*((double*)init+2*(N+1)+i);
        x[i]=*((double*)init+3*(N+1)+i);
        Ax[i]=*((double*)init+4*(N+1)+i);
        dx=*((double*)init+5*(N+1)+0);
    }

    double inital[5][N+1];
    for(int i=0;i<N+1;i++){
        inital[0][i]=u[i];
        inital[1][i]=p[i];
        inital[2][i]=rho[i];
        inital[3][i]=x[i];
        inital[4][i]=Ax[i];
    }
    double *pinital;
    pinital=&inital[0][0];

	class Laval laval;
    double *Out,*pos_UFH;
    Out = laval.macCormack(pinital, dx, CFL, Cx, errlimit, iterlimit, N);
    for(int i=0;i<N+1;i++){
        u[i]=*((double*)Out+0*(N+1)+i);
        p[i]=*((double*)Out+1*(N+1)+i);
        rho[i]=*((double*)Out+2*(N+1)+i);
    }

    double Maa[N+1],flux[N+1];
    for(i=0;i<N+1;i++){
    Maa[i] = u[i] / pow((1.4 * p[i] / rho[i]),0.5);
	flux[i] = rho[i] * u[i] * Ax[i];
    }
    std::string prefix = "./Output";
    if (_access(prefix.c_str(), 0) == -1)
        _mkdir(prefix.c_str());    
    FILE* fpout;
    char str[20];
    sprintf(str, "./Output/Laval_Cx=%.2f.dat", Cx);
    char *path=str;
    fpout = fopen(path, "w+");
    fprintf(fpout, "variables=\"x\" \"Ma\" \"Velocity\" \"Pressure\" \"Density\" \"Flux\"\nzone T=\"Laval\" i=%d F=POINT\n",N+1);
    for (int i = 0; i < N+1; i++){
        fprintf(fpout, "%15.8f\t%15.8f\t%15.8f\t%15.8f\t%15.8f\t%15.8f\n", x[i],Maa[i],u[i],p[i],rho[i],flux[i]);  //此函数把数据写入文本
    }
    fclose(fpout);

    Laval_precise(xl, xr, N, pe, Ma);
    endTime = clock();
    cout << "Calculation Finished\nDurition="<< (double)(endTime - startTime)/1000 <<"s\n"<<endl;
    system("pause");
    return 0;
}